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We construct combined electric and magnetic field variables which independently represent energy 
flows in the forward and backward directions respectively, and use these to re-formulate Maxwell's 
equations. These variables enable us to not only judge the effect and significance of backward- 
travelling field components, but also to discard them when appropriate. They thereby have the 
potential to simplify numerical simulations, leading to potential speed gains of up to 100% over 
standard FDTD or PSSD simulations. We present results for various illustrative situations, includ- 
ing an example application to second harmonic generation in periodically poled lithium niobate. 
These field variables are also used to derive both envelope equations useful for narrow-band pulse 
propagation, and a second order wave equation. Alternative definitions are also presented. 



Published as Phys. Rev. A72, 063807 (2005). A more 
detailed derivation of the fields and equations herein can 
he found at http://arxiv.org/abs/physics/0611216 



I. INTRODUCTION 

We introduce electro-magnetic field variables that 
are designed to have directional characteristics. These 
variables have the potential to speed up numerical sim- 
ulations, while providing valuable insight into the pro- 
cess of optical pulse propagation at the same time. Sim- 
ple plane-polarized versions of for a dispersionless 
medium were originally proposed by Fleck at the begin- 
ning of ref. [Q , although he did not use them in the rest 
of the paper. In the generalized form defined below, it is 
possible to use them to advantage in practical situations. 
We note that a different approach to directional pulse 
propagation based on projection operators was proposed 
by Kolesik et.al. ^ ||]; there is also the recent work of 
Ferrando et.al. [^| based on a second order wave equa- 
tion. 

The essential characteristic of G+ and G~ is that they 
represent energy fluxes directed in the forward and back- 
ward directions respectively. This implies that G+ is the 
appropriate variable to use in situations where pulses are 
travelling only in the forward direction. Indeed, as we 
will explain, optimal construction of G+ makes G~ neg- 
ligible under these circumstances, and the computational 
effort can then be halved by neglecting G~ altogether. 

If we apply a z-propagated pseudospectral spatial- 
domain (PSSD) algorithm we also gain a fast and flex- 
ible treatment of dispersion and nonlinear effects, which 
significantly outperforms standard finite difference time- 
domain (FDTD) methods 10,0]. Further, since many au- 
thors (including the recent H ) assume that the backward 
field is negligible in any case, the explicit appearance of 
G~ within our formalism provides a direct test of the 
validity of this assumption. A further advantage of G^ 
is that it is as easy to include magneto-optic effects as 
electro-optic effects such as dispersion and nonlinearity. 



It is in situations involving both electro- and magneto- 
optic effects where we achieve the greatest computational 
speed increase - potentially up to 100% faster. Moreover, 
even if one chooses to propagate an optical pulse using 
E and iJ, it is still easy to analyse its directional charac- 
teristic by constructing G^ after the event. 

After reviewing Fleck's original form of the G^ vari- 
ables at the start of section ||, we proceed to discuss 
how to represent the permittivity and permeability of 
the medium; this is a crucial step in the optimal con- 
struction of G^ in a generalized form. The treatment of 
nonlinearities and the calculation of energy and flux are 
also covered. 



In Section III, we derive a first-order wave equation, 



both in a form that is fully equivalent to Maxwell's equa- 
tions, and in a more useful one that is applicable in the 
transverse field limit. In section IV, we demonstrate a 



simple procedure for the numerical implementation of 
G^ simulations; techniques for specifying initial condi- 
tions and for handling dispersion and nonlinear effects 
are examined in detail. We take as an example the 
case of second harmonic generation in periodically-poled 
lithium niobate, to demonstrate that our method can be 
applied to practical as well as illustrative simulations. In 
all cases, we retain both G+ and G~, but show that with 
optimal construction, G~ can be made negligible. 

In section ^ we derive a propagation equation for en- 



velopes based on G^\ in section VI we develop a second 



order wave equation; and in section VII we propose al- 



ternative definitions for field variables with directional 
properties. Finally, in section VIII, we present our con- 
clusions. 



II. DEFINITIONS 

For plane-polarized fields, propagating in the z direc- 
tion in a dispersionless medium. Fleck defined the direc- 
tion field variables 



G± = ^eE^±^Hy. 



(1) 
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Their directional properties are apparent from the form 
of the Poynting vector 

S^E^Hy = -i_[G+2-G-2], (2) 



which shows that and G~ are associated with positive 
and negative energy flux respectively. Unfortunately, if 
eqn. ([^) is used to describe a forward-propagating pulse 
in a dispersive medium, it will contain significant con- 
tributions from both G+ and G^. We therefore need to 
generalize the construction in order to make the concept 
useful in practical situations. 



A. Medium Parameters 

The definitions of G* (and their generalized vector 
counterparts G^ , introduced below) depend on the prop- 
erties of the propagation medium through the permittiv- 
ity e and permeability /i. In principle it would be attrac- 
tive to define G^ using the exact values of e, (including 
the nonlinearity) , but this is usually impractical, and we 
will instead use "reference" values Er, Hr, chosen to be 
as close as practicable to the true medium properties, 
typically by including all the dispersive properties. 

In the frequency domain (indicated by tildes), we 
write 

e = iri^) + ec{uo) = al{uo) + ar{uj) aduj), (3) 

jl = jlr{uj) + ilc{uj) = PI{UJ) + fJriuj) j3c{uj), (4) 

where the correction parameters ic and jlc represent the 
discrepancy between the true values and the reference. 
The smaller these correction terms are, the better the 
match, and the more likely it is that a description involv- 
ing only G^ will suffice. Note that since the definitions of 
G^ depend on the square roots of e and /i, we introduce 
the a and /3 parameters, which will feature prominently 
(along with their time domain counterparts a,/?), in the 
generalized definitions of G^ that follow. 

By using these frequency dependent parameters in the 
generalized definitions of G^, we are able to propagate 
pulses using only the G^ variable, a gain in both math- 
ematical simplicity and computational speed. 



B. variables 

The generalized definitions of the G^ variables in the 
frequency and time domains are 

G^{uj) = ar{Lj)E{uj)±ux j3r{uj)H{uj), (5) 

G°{u) = u- [/3r(w)i?(w)j ; (6) 

or G^{t) = ar{t)*E{t)±ux l3r{t)*H{t),{l) 

G°{t) = u- \pr{t)*H{t) \ , (8) 



where u is the unit vector in the direction of propagation, 
and ar{t) and I3r{t) are the (inverse) Fourier transformed 
versions of ar{uj) and (3r{i^)- The symbol is used to 
denote a convolution: a*b = J a{t — t')b{t)dt' . The vari- 
able G° involves the longitudinal part of the magnetic 
field, which is eliminated in the u x H operation in eqn. 
(^. Although we will generally make a transverse ap- 
proximation in which G° = and u ■ G^ = 0, we retain 
the longitudinal parts of the field to ensure a complete 
description. To avoid cluttering the notation, we do not 
apply tildes to the spectral forms of the field quantities 
G^, G°, E, H, and rely instead on the arguments {t or 
u) or the context, to distinguish between domains. 

Inverting eqn. (^ gives the following expressions for 
the electric and magnetic fields as a function of G^ and 
G° 



E{lu) 



1 



1 



G+ {uj) + G- {uj) 
G+{uj)-G-{u) 



uG°{lo) 



(9) 



(10) 



The divergence of G^ , allowing for both charge density 
p and current density J, is 

V.G± = %p±Z^PrU-{G+ + G-)TPrU-J{i^) 
a'^ 2 Ur V / 

We note that this is zero when the p and J are zero, as 
long as there is no longitudinal electric field. 

Different choices of er,P-r produce different G^ pairs. 
Whilst using the true values to describe a forward prop- 
agating pulse results in G^ = 0, any other choice of 
reference will produce a non-zero G~ component that 
co-propagates with G"*". Note that this G^ still has an 
energy flux directed in the reverse direction (— u), but 
travels forwards with the G^ with which it is tightly cou- 
pled. 

We will almost always choose , Jlr to include the en- 
tire linear dispersion of the medium. We exclude the 
nonlinearity because it removes the ability to reconstruct 
-B, H fields uniquely from the G^ , as can be seen from 
eqns. (^, ^, ^ |l^), which will become nonlinear in E and 
H. 

The vectorized definitions of G* accommodate any po- 
larization of the E and H fields. For propagation along 
the z axis, the x component of G^ (G^) will contain E^ 
and Hy] and similarly G^, will contain Ey and H^. It 
is then a simple matter to see how linearly or circularly 
polarized E and H fields can be represented in terms of 
G^ . The definitions are also easily generalized to include 
birefringent media, provided the propagation direction 
and transverse coordinate axes are such that and p,r 
become diagonal matrices. 

Finally, note that G^ bear some resemblance to Bel- 
trami variables (see e.g. ^, |lOj) which are defined 
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as Q = \ftE + i^fJlH] but they differ in two important 
respects. First, a given Beltrami Q defines E and H 
uniquely, whereas both G"'" and G~ are needed to do the 
same. Secondly, Q does not assume any preferred direc- 
tion, whereas the variables include the direction u in 
their definition. Further, Beltrami variables are not de- 
fined using the full time (or frequency) dependence of e, /i 
as we use for G^ in eqns. (^, although presumably 
this would be possible. 



C. Nonlinearities 

Since it is usually impractical to include nonlinearities 
in the reference parameters, these will normally appear 
in the correction terms ec, ^c- As an example, consider 
a n-th order (electric) nonlinearity, in which case €c{t) — 
x(")(i) and 



a,{Lo) = [5,(c^)]-\jM")(t)*i?(t)"-i 



(12) 



a,(0 = :r-i{[5,(u;)]-\x^")(a;).:r[£;(t)"-i]}a3) 

where 3^[...] is the Fourier transform (FT) from time to 
frequency, and E(t) can be found from eqn. (|^). If the 
reference parameters Ur contain dispersion (which will 
be the typical case), we can see from eqn. ( p^ ) that 
this will make adoj) dispersive even if x^"^ is instanta- 
neous. In the case of an instantaneous nonlinearity, this 
adds more computational work (an extra two FTs), al- 
though for non-instantaneous ones we needed the FTs 
anyway. If the nonlinearity is instantaneous and the 
reference parameters are non-dispersive, we have simply 
af ^(f) = a-".x(").2-"+i [G+ + G"]""'. 



D. Energy and Flux 

The G* are intrinsically directional and do not rely 
on a carrier wave to impart their directionality. This 
becomes clear when the Foynting vector is expressed in 
terms of G^. For transverse fields and dispersive refer- 
ence parameters, we obtain 



interpretation that for particular E and H fields, G"^ cor- 
responds to the energy flux directed forward (along u), 
and G^ to flux directed backward. The need for this 
distinction between the direction of the flux due to a G* 



field, and its direction of travel has already arisen in [I B 
above. 

We can also calculate the energy density of the EM 
field, U{t) = \e * E{t) ■ E{t) -h i/it * H{t) ■ H{t). For 
transverse fields and a non-dispersive reference, while still 
allowing for medium dispersion, the energy density in 
terms of G^ is 



U 



Mr 



*G^ 



± _ ii 

e,. Mr 



•G^ + i 
* G+ 1 • G" 



Mr 



* G- 



G+. 



G- 



(17) 



S 
S 



E X H 
1 



- (j-i * G-) • {j-' [p; 



* G^ 



G- 



(14) 



u. (15) 



Notice the cross terms, which appear whenever there 
is a mismatch between the reference and medium param- 
eters. These occur because of the interference between 
the G+ and G~ contributions to the field. 

For a dispersive reference, the relevant formulae for S 
and U are relatively complicated because of the appear- 
ance of cross terms and/or convolutions. However, this 
should not produce a significant overhead in numerical 
simulations because the code will be switching between 
time and frequency domains at each step, allowing S and 
11 to be calculated in whatever way is most efficient. 



E. Co-moving frame 

We now consider using a moving reference frame. This 
is particularly useful in a space-propagated model where 
the pulse is held as a function of time, since it will stay 
nearly centered when propagating forwards. A simple 
choice of frame speed might be the phase velocity at the 
centre frequency of the pulse, which minimises the mo- 
tion of the carrier-like oscillations; however, the pulse as 
a whole will move within the frame because of its dif- 
ferent group velocity. The frame translation for a speed 
Cf = 1/afPf is 



(18) 
(19) 



For dispersionless reference parameters, this becomes 
simply 



where 7 is the distance travelled in the direction of u. 
Thus 



1 



G+ ■G+ -G- ■ G- 



(16) 



dt = dt 
V 



V - -dt. 



(20) 
(21) 



Since both the G^ • G^ terms are real and positive, 
we see that G^ and G~ contribute positive and nega- 
tive energy fluxes respectively. This leads to the simple 



In vector calculations, we need to know how this frame 
translation transforms the curl and divergence opera- 
tions. The divergence is a straightforward consequence 
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of eqn. ( |21| ) , and the curl of an arbitrary vector Q trans- 
forms to 



VxQ = V'xQ - afPfU x dtQ. 



(22) 



The ratio of the reference speed (the phase velocity in 
the reference "medium" described by e^j/ir) and frame 
speeds is 



(23) 



If we choose to give aj and (3f a frequency depen- 
dence, we have defined a "dispersive frame" , where dif- 
ferent frequency components travel at different speeds. 
In such a frame, any matching dispersive evolution (i.e. 
where af = ar and f3f = (3r) results in no change to the 
pulse profile. However, at the end of the simulation, we 
need to transform from the dispersive frame back into a 
normal (non-dispersive) laboratory frame. Moreover, us- 
ing a dispersive frame can give rise to numerical stability 
problems. 



III. FIRST ORDER WAVE EQUATION 

We now derive a set of first-order differential equa- 
tions for the forward and backward directed fields G^, 
and use the moving frame set out above in eqns. (2l|, 
p^ ). We assume that the medium is continuous, so that 
dz^ = dz^i = 0, where dq = d/dq. This does not impose 
a significant restriction in practice, since a simulation 
propagated forwards in space can easily handle interfaces 
between different media. 



A. Derivation 

For a vector derivation of propagation equations 
for G^, we start with the two relevant (source free) 
Maxwell's equations. Writing them in frequency space, 
with (5^ = e and — p,; and taking the cross product of 
u and the \7 x H equation yields 



u X 



V X iJ 

V X £: 



-lujdPu X E, 



(24) 
(25) 



We now separate the correction components (depend- 
ing on etc, Pc) from the reference components (depending 
on ctr, Pr), and substitute expressions containing G^, G° 
by referring to eqns. (||) and (10). We also note that the 
terms involving G^ decouple from those involving G° . 
Hence 



V X G — ^lUOirPr U X G 
lUJOLcBr ^ 

T — - — u X 



2 

ILVOtrPc - 



-U X 



G+ + G- 



G+ -G- 



VG° — ±lUjar(3rUG° ± lUJCtrPcUG" 



(28) 
(29) 



For media whose magnetic behaviour is matched per- 
fectly by the reference parameters (i.e. /3c = 0), this 
simplifies to 



V X G^ = ^lUJarPr U X G^ =F 
VG° = ±lUjari3rUG° . 



G+ + G- 



(30) 
(31) 



For propagation along the z axis, in the plane polarized 
(_Ej., Hy) limit, the curl becomes dz, and G^ is replaced 
with G^. In the transverse field case, eqn. (|l|) (or (|2g|)) 
can be ignored. 

The equations above are written to suggest a spatially 
directed propagation (along u), and indeed are most 
straightforwardly solved that way. However, a simple 
rearrangement of the terms leads to a i-directed propa- 
gation model, although, in its present form, the time-like 
evolution of eqn. (^ is obscured. Whilst t-directed 
propagation has some advantages, it makes the treat- 
ment of dispersion and other time-memory effects more 
demanding, as discussed by Kolesik and Moloney |^ and 
TyrreUet.al. @. 

In the time domain, eqn. (pO|) becomes 



V X G± ±dt 
±dt 



{ar * Pr) * (^U X G^^ 

ac * p., 



G+ + G- 



(32) 



Multiplying respectively by Pr and otr and taking sums 
and differences leads to 

V X UrE ± ux {v X PrH^ = +iujarP^ H 

T lujPrCt^u X E. (26) 

Noting the similarities between this and eqn. (|^), we 
now reorganize using standard vector identities for V x 
{A X B) and u x {u x A). Finally, we arrive at a curl 
equation for G^, namely 



Note that reversing the direction of propagation by 
changing u to —u reverses the roles of G+ and G^. 



B. Longitudinal E 



V X G^ 



^lUJ I PrCi^ U X E ± arj?U X 

+ iuar~p^uG° tvg°. 



u X H 



In the same way that the construction of G^ ignores 
the contribution from H along t he p ropagation vector u, 
and forces us to define G° , eqn. (gj) ignores information 
about the time-evolution of the longitudinal part of E. 
We can rectify this by using V • (u x _ff ) = m • (V x _ff), 
to get 



(27) 



Q!,.V • G^ — G — —luia'^PrU 



G+ + G- 



(33) 
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This is the difference of the source-free divergences for 
calculated in eqn. (|l|). Thus eqns. (||, |2|, 
|33| ) provide another way of solving the complete set of 
source-free Maxwell's equations using an alternative ba- 
sis. Clearly, however, our basis of G^, G° is most useful 
for fields propagating mainly along one axis (i.e. u), par- 
ticularly in the limit of transverse fields, where only eqn. 
(^) needs to be solved. 



C. Co-moving frame 



We can transform eqn. ( |28D directly into a moving 
frame using eqn. which gives 



V X G± 



Ti^^ar(3r (1 T ^ ^ 

G+ + G- 



T — - — u X 



2 



G+ -G- 



(34) 



Here we have not shown the transformed (non- 
transverse) eqns. (^ p3| ) in the interest of brevity, but 
they are easy to calculate if needed. 

One nice property of this equation is that matching the 
frame velocity to the phase velocity causes the carrier-like 
oscillations in the forward travelling fields G"*" to freeze 
in place, leaving only the evolution due to the correction 
terms. If we are prepared to make the common assump- 
tion of only forward-travelling pulses, we will have man- 
aged to greatly reduce the rate of change of the fields. 
This in turn will allow coarser numerical resolutions to 
be employed in numerical simulations, leading to signifi- 
cant speed advantages over and above those obtained by 
assuming G~ = 0. 



E. Decoupled Wave Equations 

We can make the most of our approach by decou- 
pling G^ from G^, enabling the two first-order coupled 
Maxwell's eqns. ( p^ , p5| ) or G^ eqns. ( |2^ ) (or co-moving 
form eqn. (|3J)) to be reduced to two uncoupled first- 
order equations. The equation describing propagation in 
the "uninteresting "direction can then be discarded, leav- 
ing one first-order equation where there were originally 
two. 

This step requires an approximation, although since we 
can perfectly match the reference parameters (a,-, (3r) to 
the material dispersion, it is not a very stringent one. 
Since the correction parameters etc, (3c depend only on 
nonlinear effects, they will in general be small, keeping 
cross coupling between G'^ and G~ minimal. Further, 
whilst the G+ field will rotate forwards according to its 
wavevector (i.e. e"''*'^^, with k = ar(3rUj), the G^ field will 
rotate backwards at the same rate (i.e. at 6"''=^). This 
means the correction terms for G"*" will contain both an 
in-sync component from G"'", and a component from G^ 
with a large detuning. Since this detuning (amounting 
to e"^**^^) will usually be large compared to the spatial 
bandwidth of the pulse, we can apply a rotating wave 
approximation and average the G~ contribution to zero. 
After applying the same steps to the G~ equation as well, 
eqn. (|2q) becomes 



V X G± 



^lLL}ar(3r U X G^ 

T — ^ — u X G± T — — u x G±,(35) 



IV. SIMULATING 



D. Time vs space propagation 

In FDTD solutions of Maxwell's equations, optical 
pulses travel either forwards or backwards in space as 
they propagate (or march) forward with time. However, 
most nonlinear optical simulations are done in a space- 
propagated picture; with the consequence that optical 
pulses travel either forwards or backwards in time as 
the calculations propagate (march) through space. 

Since we follow the space-propagated picture, the pulse 
travelling forward in time will be described by G+, and 
the one travelling backward by G~ . 

Note that any backward travelling pulse in a z- 
propagated picture is travelling backwards in time while 
propagating forwards in space. Although at first this 
might seem non-causal, it is in fact the way that the 
simulation represents a pulse which we would normally 
describe as propagating backwards (i.e. in the direction 
—u). This is clear from the wave equations; swapping the 
sign of the propagation direction u swaps the behaviour 
of G+ and G". 



We now examine the procedure needed to simulate 
wave propagation using the G^ variables. This will clar- 
ify various practical issues as well as illuminate some of 
the less-obvious features of our approach. We consider a 
plane-polarized EM wave propagating along z in a non- 
magnetic medium with dispersion and a weak nonlinear- 
ity. Since our aim is to explain the fundamental principles 
of the use of G^ variables, we first present a number of 
simple examples. 

Our numerical simulations of the wave equations 
are implemented by straightforward adaption of the 
PSSD technique @. In PSSD, fields are stored as func- 
tions of time, and fast Fourier transforms (FFTs) are 
used to convert to the frequency domain for the calcula- 
tion of pseudospectral derivatives and the effects of dis- 
persion. This technique allows the simple application of 
arbitrary dispersion, which becomes a simple multipli- 
cation in frequency space. Fields are then transformed 
back to the time domain, where the nonlinear effects 
are calculated, before propagating the fields forward in 
space. Computational details, such as how to design the 
mesh and control the accuracy of the simulations are well 
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known (e.g. the Courant and Nyquist criteria), and can 
be found in a range of sources (e.g. 13| 1^). 

When apphed to fields, the basic spatially- 

propagated PSSD algorithm does not change, but the 
wave equation to be solved is now eqn. (^) instead of 
Maxwell's equations. This means that the full flexibility 
of PSSD is harnessed with the advantages of fields to 
give a powerful and efficient combination. 



A. Simulation speed 

The computational speed of any PSSD-type propaga- 
tion depends primarily on the time spent doing FTs. In 
the PSSD technique described in ||], five FTs are used, 
two forward and back pairs, and one (forward only) to 
calculate the derivative of the electric displacement D. 
If magnetic dispersion were present, PSSD would require 
an extra FT for the magnetic induction B, making six 
FTs in all. 

In contrast, a simulation requires only three FTs. 
This comprises two forward FTs which are used to cal- 
culate the derivative for the dispersion and nonlinearity, 
and one backward FT is used to change back into the 
time domain; the two derivatives are combined in the fre- 
quency domain where the problem becomes linear. Such 
a simulation will therefore run 67% faster than the cor- 
responding E and H PSSD algorithm; or 100% if there 
is also magnetic dispersion. 

To include both and would require six FTs; one 
more than the usual PSSD case, but the same if magnetic 
dispersion needs to be included. 



B. Implementation 

As a first step, we divide the total permittivity into 
three: a reference component with constant permittivity 
ir, a linear dispersion correction e^(w), and an instanta- 
neous nonlinearity e^^; the permeability has the vacuum 
value /iQ. The medium properties can therefore be rep- 
resented in the following fashion 

e» = e~, + ?f(cc>) + gf^ (36) 
= al + &r&c i^) +(^ra^^ (37) 

This particular breakdown of e is for illustrative pur- 
poses only; in practice, we would choose a dispersive ref- 
erence and try to leave only nonlinear terms in the cor- 
rection parameters (i.e. use e(w) = er{u>) + e^^)- We 
might also regroup various terms to optimise the numer- 
ical performance. 

The first order evolution equation, specialized from 
eqn. (|30|), is 



^Pr (1 T C) Gt T [Gt + G-] 



[Gt + Gx ] 



(38) 



where the RHS contains respectively a reference carrier 
term (cx cXr), a linear dispersion term (oc a^), and a 
nonlinear polarization term (cx a^^). We integrate for- 
ward in z using a split-step method, where each term is 
integrated through dz in sequence. This procedure is ac- 
curate to first order, so we need to ensure Sz is sufficiently 
small. 

In this simple case, the reference term merely applies 
a complex rotation to the field in frequency space repre- 
sented by 



Gxi — G't{z)xeyip ^lujarPr (l T -^^ 



(39) 



If the frame velocity is chosen to be the same as the 
phase velocity given by the reference parameters, the 
Gt field no longer undergoes any reference evolution as 
it propagates. This contrasts with the usual approach, 
which is to match the frame velocity to the group veloc- 
ity. However, with our choice of reference parameters, the 
group velocity corrections appear in the second RHS term 
as part of . If we were to propagate in a group ve- 
locity frame, we would retain part of the reference term, 
which would then cancel with part of the group velocity 
contribution from the dispersion term. This could lead to 
a better overall cancellation, just as in the usual E field 
approaches. If that were our aim, we could indeed easily 
rearrange eqn. ( p8| ) to incorporate such a cancellation, 
and then solve the equation appropriately. 

The next step is to solve for the linear dispersion — 
far- Fortunately, this part of the equation is also easy 
to solve exactly in the frequency domain, through the 
operation 



Gx2 = X exp [Ttk G^^^.Sz 
T 7, \ — V-'xl 



G ^ 



(40) 



Although both reference and dispersion steps can be 
solved using exponentials, there is an important differ- 
ence. The reference evolution of G"*" depends only on 
G^ , whereas the dispersion evolution depends on the sum 
G+-f G~, since the dispersion acts on the electric field. In 
a forward-only approximation where G~ = 0, it is trivial 
to combine these first two steps, as in most approaches 
to solving for the propagation of optical pulses. 

The third and final step is performed by transforming 
into the time domain and solving for the n-th order non- 
linear effects. Since the reference ar,(3r are constants. 



-^MjTjn a simple Euler method gives 



G^iz + Sz) 



G: 



An) 



x2 



± 



2" 



'.dt [Gi 



G ^ 



-1 
oz. 



(41) 



For a narrow-band field centred at c^o , the time deriva- 
tive would be dominated by (and proportional to) wq. In 
most envelope theories we see only this factor luq in the 
analogous expression; although correction terms exist for 
wider-band fields 0, 
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Initial conditions: matching a pulse to the 
medium 



Most descriptions of pulse propagation start with ini- 
tial conditions chosen to represent a pulse travelling for- 
ward in the medium. Here we consider how to choose the 
best initial conditions for in the case of a pulse trav- 
elling only in the forward u direction. They are based on 
the best practical parameterization of the medium ei{uj), 
fli{uj), which need not be the same as and fir- Assum- 
ing only the electric field E{llj) of the pulse is known, the 
procedure is: 

(1) Choose ii and fli to be as close as possible to the 
actual medium parameters e, jl. One might even try to 
put the nonlinear properties into ii and fli as well, but 
only if one can get a solution for steps (2) and (3) below 
with this added complication. 

(2) Calculate H{uj) corresponding to E{uj) for a forward 
traveUing pulse, so that a G~ based on ii, fli would be 
zero: 



Eiu). 



(42) 



(3) Calculate an initial using the chosen reference 
parameters er{Lo) and i.ir{Lo), given our initial E{lj) and 
H{uo) fields: 



E{lo). (43) 



Note that step (3) is unnecessary if — and jli — jlr- 
Notwithstanding step (2), G~ is only eliminated from 
the simulation if both q and jli are perfect matches to 
the material parameters. If the {li, jli) values are good, 
but {ir, jlr) less SO, a weak G~ field will co-propagate 
forwards with G^, even though the Poynting vector of 
the G' is directed backwards. If (e^, jii) is a bad match 
as well, the initial G~ will have a component that travels 
backwards, its magnitude corresponding to that of the re- 
flection between a medium with parameters (e^, jli) and 
one with the actual parameters. Since it is usually pos- 
sible to include all the linear dispersive properties in e^, 
jli, any discrepancy is likely to be due to the nonlinear 
contribution, and consequently very small. 

Figure |l| shows how different choices of reference pa- 
rameter affect the G^ fields required to model a simple 
forward-propagating few-cycle pulse at 500nm in fused 
silica. Note that although the G^ fields in (a) and (b) 
are directly proportional to E and in (c), the use of 
a dispersive reference means that a deconvolution would 
be needed (if in the time-domain) to transform from G* 
to E and H. 

In figure [^(a), the mismatch between the reference pa- 
rameters (with n = 1) and the actual medium (n « 1.5 
at the 500 nm pulse center wavelength) causes a signif- 
icant co-propagating G~ component to appear; this is 
improved in (b) where the reference parameters specify a 



constant refractive index close to that at the centre fre- 
quency of the initial pulse. Since the mismatch between 
the reference and the true material properties is due only 
to the material dispersion, the initial co-propagating G~ 
component is smaller in figure ^(b) than in figure 0(a). 
The reduction in the size of G~ is rather smaller than 
might be expected, mainly because although fused silica 
has a refractive index of about 1.5 at 500nm, we have 
used non-dispersive reference parameters with a refrac- 
tive index of 1.5 at all frequencies. 

In figure |l(b), although the construction of G* has 
the phase velocity reasonably well matched, the group 
velocity of the pulse is poorly matched. In addition there 
is a smaller effect caused by the wide-band nature of the 
pulse, where the reference parameters are (even) less well 
matched to frequency components away from the center 
frequency. 
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FIG. 1: A 500nm pulse in fused silica, represented with (a) a 
vacuum reference, (b) a fixed refractive index reference (c) a 
perfectly matched dispersive reference. In all cases the initial- 
ization parameters are those that perfectly match the disper- 
sive properties of the medium. Solid line: field. Dashed 
line: G~ field. 

The conclusion is that, for any pulse propagating in 
a material whose dispersion is not perfectly matched 
by the reference over the pulse bandwidth, a finite co- 
propagating G~ will appear. This will be made up of fre- 
quency components whose phase velocity in the medium 
do not match the phase velocity given by the reference. 
Thus, in typical dispersive media, only very narrow-band 
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pulses result in a negligible G~ for non-dispersive refer- 
ence parameters. However, the mismatch between the 
reference and the true material properties can be com- 
pletely removed by using a dispersive reference identical 
to that of the material being simulated. The results of 
this are shown in figure |l|(c) , where G~ is identically zero. 



D. Dispersive propagation 

We now present a variety of numerical results demon- 
strating pulse propagation in a dispersive medium. We 
take the medium to have the properties of fused silica, 
but we do not include nonlinear effects for the moment. 
Our aim is to give a flavour of what the fields look 
like for different reference parameters. The choice of ref- 
erence is important because, as explained earlier, if the 
reference is not perfectly matched to the actual medium, 
a forward travelling pulse will contain a G~ wave co- 
propagating with the main component. Usually we 
will want to choose a reference that makes G~ negligible, 
so we can save computational effort. 

We deliberately choose ultra-short pulses containing 
only a few optical cycles to demonstrate the flexibility of 
our method in the short pulse limit. 

Figure ^ shows the results for the fields in fig. |^ after 
propagating 15/im in fused silica with the nonlinearity 
ignored. In all cases, the initial size of G~ is broadly 
maintained and, in particular, it remains zero when the 
reference parameters are perfectly matched. Although 
the fields in ^(a) and ||(b) are directly proportional 
to E and H, in (c), the use of a dispersive reference means 
that in that case a deconvolution is needed to transform 
from G± to E and H. 

We can also consider the effect of neglecting a finite 
(but significant) G~ field, where the G+ part of the pulse 
then undergoes the wrong dispersion. This is because the 
dispersive correction part (see e.g. eqn. (|3|)) depends on 
G+ 4- G^ , and thus, without G~ , will be either too big or 
too small. This problem is avoided by using a dispersive 
reference identical to that of the material being simu- 
lated. The results of this are shown in figure ||(c), where 
G~ is always identically zero and no approximation is 
necessary to omit G~ . 

In figure |^ we show the result for a simulation with 
both a perfectly matched dispersive reference and a per- 
fectly matched dispersive frame. Since all the material 
properties are included in the reference parameters, and 
we pick a frame that exactly matches the propagation, 
fig. 1^ looks identical to the initial state in fig. ^(c). We 
can recover the expected lab-frame final state by trans- 
forming fig. js] out of its dispersive frame, and so get a 
graph identical to fig. ||(c). 

The main message from these simulations is that the 
better matched the reference parameters are to the ma- 
terial parameters, the smaller the co-propagating G^. 
For a perfectly matched reference, the co-propagating 
G~ vanishes. Also, the better matched the frame is to 
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FIG. 2: The same pulse as in figure |l|, after propagating 
15/xm; represented in (a) a vacuum reference, (b) a fixed re- 
fractive index reference (c) a perfectly matched dispersive ref- 
erence. Solid line: G+ field. Dashed line: G" field. 
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FIG. 3: Simulation results showing G for perfect reference 
(tr ~ ^3aica{i^)) and a matched dispersive frame. 



the material parameters, the slower the evolution of the 
pulse shape. However, we then have to do more work 
to transform the final state of the pulse (in its moving 
frame) into the stationary-frame counterpart we would 
see in the lab - although for a linearly dispersive frame, 
the transformation is straightforward. 



E. Nonlinear propagation 

We now demonstrate some simple pulse propagations 
in nonlinear media. Since neither the initial conditions 
(determined by nor the reference parameters (de- 
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termined by er,fJ.r) include the nonlinearity, the pulse 
is not perfectly forward propagating, and a small "re- 
flection" occurs as the pulse starts propagating in the 
nonlinear medium. 

Figure ^ shows how pulses similar to those in fig. |l| 
look after propagating 10/im through fused silica. The 
pulse parameters were adjusted to give a clearer final 
pulse shape. We see the same pattern as in figs ^ and ||, 
where a weak G~ remains except for perfectly matched 
reference parameters. Note, however, that the addition 
of nonlinearity does not cause the size of the G~ field to 
change significantly during propagation. 
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FIG. 4: A similar pulse as above, after propagating lO/im 
through fused silica; represented in (a) a vacuum reference, 
(b) a fixed refractive index reference {n = 1.5), (c) a per- 
fectly matched dispersive reference. Parameters have been 
adjusted to give a clearer final pulse shape. Solid line: G""" 
field. Dashed line: G~ field. 



We now apply our approach to the practical problem 
of second harmonic generation in 120/xm of periodically 
poled lithium niobate. The results are shown on figure 
H, and agree with the simulations of Tyrrell et.al. Q. 
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FIG. 5: Second harmonic generation in 120 fim of LiNOa, 
periodically poled at 6.05/im. Clockwise from top left: initial 
pulse, final pulse, second harmonic power, final pulse spec- 
trum. Solid line: G"*" field. Dot-dashed line: G~ field. 



F. Some remarks on layered media 

A complication arises when propagating a pulse 
through layers of material with significantly different dis- 
persions. Because the definitions are carefully con- 
structed to match the propagation medium, variables 
ideal for one layer (and so ensuring G~ = 0) will not be 
ideal for another. This gives us two options: (a) either 
retain the G~ field in the description, or (b) at each layer 
boundary, switch to a set of G^ variables matched to that 
medium. Option (a) is simpler, but it is not necessarily 
computationally efficient and leads to complications in- 
volving refiections from the interfaces. Option (b) is more 
efficient computationally when we are only interested in 
the forward-going pulse, as the effort involved in switch- 
ing G^ definitions is comparable to only a single spatial 
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step in the ongoing propagation calculation. 



G. Justifying the forward-only approximation 

The ability to accurately incorporate dispersion into 
our reference permittivity allows great control over the 
magnitude of the G~ field. Our tests have shown that 
we can confidently neglect G~ if our construction of 
accurately includes the medium dispersion, although pos- 
sible exceptions may occur in cases involving extremely 
strong nonlinearities. 

This can be seen in the case of periodically poled 
lithium niobate discussed above (see fig. ||), where the 
ratio of the G^ to G^ intensities was 1 : 10^. An even 
more rigorous test of G+'s ability to accurately simulate 
short pulse propagation was our recent study of the ef- 
fects of dispersion on carrier shocking jl^. Despite the 
strong nonlinear effects, and significant distortion to the 
pulse profiles, G^ simulations consistently produced re- 
sults in agreement with PSSD - whilst still only requiring 
half the computational effort. 

The ability to accurately model pulse propagation us- 
ing only G"*" after carefully choosing a reference permit- 
tivity clearly justifies neglecting G~, which in turn sim- 
plifies numerical simulations. 



V. ENVELOPE PROPAGATION EQUATION 

When computing the interaction of narrow-band fields, 
it is common to remove chosen carrier frequencies, and to 
evolve the envelopes rather than the complete EM fields. 
In fact, if sufficient care is taken with the approxima- 
tions, and the system simulated is well behaved, even 
quite wide-band pulses can be successfully modelled in 
this way. 

We can use an envelope approach with the G^ vari- 
ables. However, a full model requires four envelopes to 
describe the G^, just as in a complete Maxwell theory 
where envelopes are needed for both the backward and 
forward travelling E and H. A full expansion of G^ into 
forward and backward envelopes G^ , G^ would be 



waves are now eliminated, we can propagate pulses effi- 
ciently in a moving frame. This is important, because 
the backward parts in a moving frame move at twice the 
frame speed. In a full (non-envelope) simulation, we need 
somehow to filter out the backward components, as oth- 
erwise the hoped-for numerical gains are lost by the fact 
that a finer z-step is required for accurate integration. 

The first order wave equation for the forward-travelling 
envelopes defined above is 



d,5 



f _ T^(^ci./3.-fco)s^T^^{s^ + SJ*}(.4£ 



Here Gc + a* — a^, which is simple in the case of dis- 
persion but, in the presence of nonlinearity, will be the 
appropriately carrier-matched, positive frequency part of 
the permittivity correction parameter. 
If ujQ = koCr, we have 

a,g± = ^,5^^,,(^_^o)g±^!^^|g± + gJ *|(.46) 

In a suitable narrow-band limit, we should be able to 
ignore the first term on the RHS of this equation, leaving 
the evolution of the envelopes to be controlled solely by 
the correction term. The description can be easily gen- 
eralized to cases involving multiple components centred 
on different carrier frequencies. Note that this is a first- 
order envelope equation, and, as such, does not require 
the various extra approximations needed when deriving 
an envelope propagation equation from the standard (E 
field) second order wave equation. 



VI. SECOND ORDER WAVE EQUATION 



+g±(cc> T cc>o)e±*^«^ + *(co T c^o)e^''=°^ (44) 

where we have suppressed the z argument on the enve- 
lope functions for brevity. Note that the forward-like 
G~ contribution (i.e. Sj) needs a backward-travelling 
carrier, as otherwise it is not possible to match the refer- 
ence evolution terms for both Sjr and SJ ■ When inserted 
into the wave equations, this expansion results in a large 
number of terms, even for the relatively simple case of 
a third-order nonlinearity. However, we can specialize 
to the case where only forward-travelling waves are con- 
sidered, and set = 0. Since the backward-travelling 



In section HI we derived first-order wave equations for 
the field variables G^. However, since many pulse prop- 
agation theories start from a second-order form, we have 
also derived a second-order propagation equation. We 
apply the usual restriction to transverse-only fields, and 
split the medium properties (i.e. the permittivity and 
permeability) into a reference part (with Cr = 1/arPr), a 
linear dispersive part (controlled by q;^,/9^) and a non- 
linear electro-optic polarization part (P = a^a^^ * E). 
The time-domain wave equation for a non-dispersive ref- 
erence is 



2r<± 



V^G 



1 



dfG 



--dt I —dt T u X Vx 

2 ^\cr 



G+ + G- 

_ 1 
2ar 



G+ -G- 



dt 



— dt =F u X Vx 

Cr 



This is similar to the usual second-order equation for 
the electric field, but has the addition of a curl operator 
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applied to the dispersion and polarization terms. This 
second order wave equation can be solved with the use 
of an envelope-carrier representation for , as is often 
done with the standard equation for the electric field E. 
Such a derivation can be found in JTot , which contains 
both SVEA and GFEA || versions for both E and G^. 
The most general form of eqn. (E^) appears in |17|. 



±VF° = +lUjarPrU F° + lUacPrU F° {h2) 

V-{f+-F-^ = -iLoar {f3r + Pc) u ■ (^F+ + p-'j 

+ iPr + f3c)u- J. (53) 

VIII. CONCLUSIONS 



VII. ALTERNATIVE DEFINITIONS 

Just as one may decide to propagate the D field instead 
of the E field, so directional field variables in the style 
of G^ can be defined in a number of ways. Continuing 
with the pattern of combining transverse field compo- 
nents with a cross product, alternative directional fields 
are 



G'± = M X arE + j3rH, 



F 



G'° ^u- arE; (48) 
F° ^u-f3;:^B; (49) 
F'° ^u-a-^D.{50) 



The G^ or G'^ variables will best suit problems de- 
fined in terms of E and H; the are best suited to 
electric media, and the G'^ to magnetic media. In con- 
trast, the F^ or F'^ variables are more suited to D and 
B. All these definitions can be used to generate wave 
equations, by a similar procedure to that in section [11. 
A point to note is that if the wave equations are general- 
ized to include source terms, the G^ and G'^ forms (or 
F^ and F'^ forms) of the wave equations look somewhat 
different. 

As an example, here are the full first order wave equa- 
tions for the F^,F° form, which is conceptually closest 
to the UPPE (unidirectional pulse propagation equation) 
of Kolesik et.al. B, H] based on projections of I? - 



V X F"^ = ^lujarPr uy. F^ 

ILOacPr ^ 



F+ +F~ 
F+ -F 

± UX {(3r+ Pc) J, 



(51) 



We have introduced generalized forms of the direc- 
tional field variables first envisag ed by Fleck||]. We 
have demonstrated that they are associated with en- 
ergy fluxes in the forward and backward directions. 
They provide the ideal basis for the standard "forward- 
only" pulse propagation model, both improving our in- 
sight into pulse propagation, and allowing the backward- 
propagating component to be efficiently discarded if de- 
sired. By developing the theory in frequency space, we 
have shown how the dispersive properties of the propa- 
gation medium can be incorporated. 

We have derived first-order wave equations for G^ 
that are equivalent to Maxwell's equations. If disper- 
sion is included carefully, the equations decouple and 
we can get a single equation for forward-only propaga- 
tion, and hence achieve significant speed gains over di- 
rect Maxwell's equation solvers for E and H (such as 
PSSDH or FDTDg]). We have also presented a number 
of simulations demonstrating their use. 

Since the G^ variables are not restricted to use in first 
order wave equations, we have also presented an envelope 
theory and a second-order wave equation analogous to 
those regularly used in pulse propagation work. Either 
of these equations can be used to extend the practical 
applications of G^ variables into the long-pulse narrow- 
band regimes. Further, can still be constructed from 
the E and H field obtained in traditional simulations, 
enabling their use for either diagnosis or analysis. 
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